PxE

Systematic Associations Between Phenotypes and Environmental Exposures: Benchmarking Exposomic Research

NHANES

if(exists("con")) {
  dbDisconnect(con)
  remove(list=ls())
}

library(tidyverse)

## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.3.6     ✔ purrr   0.3.4
## ✔ tibble  3.1.8     ✔ dplyr   1.0.9
## ✔ tidyr   1.2.0     ✔ stringr 1.4.0
## ✔ readr   2.1.2     ✔ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()

library(DBI)
library(ggsci)
library(DT)
library(ggrepel)
library(cowplot)
library(reactable)
library(gt)
library(broom)

ep_assoc_summary_across_models <- function(summary_stats, glanced_stats) {
  summary_stats_wide <- summary_stats |> pivot_wider(names_from = "model_type", values_from = c("estimate", "std.error", "statistic", "p.value")) 
  summary_stats_wide <- summary_stats_wide |> mutate(estimate_diff = estimate_adjusted-estimate_unadjusted)
  adj_vs_base <- glanced |> select(-c(adj.r2, df.residual, null.deviance, df.null, deviance)) |> pivot_wider(names_from=model_type, values_from = c("rsq", "nobs", "AIC", "BIC"))
  adj_vs_base <- adj_vs_base |> mutate(rsq_adjusted_base_diff=rsq_adjusted-rsq_base, rsq_adjusted_diff = rsq_adjusted-rsq_unadjusted)
  summary_stats_wide |> left_join(adj_vs_base, by=c("evarname", "pvarname", "exposure_table_name", "phenotype_table_name"))
}

remove_units_from_string <- function(vardesc) {
  gsub("\\(.*\\)$","", vardesc)
}

con <- DBI::dbConnect(RSQLite::SQLite(), dbname='../db/pe_summary_stats.sqlite')
varnames <- tbl(con, "variable_names_epcf")
adjusted_meta <- tbl(con, "adjusted_meta")
unadjusted_meta <- tbl(con, "unadjusted_meta")
adjusted_meta <- adjusted_meta |> left_join(unadjusted_meta |> select(evarname, pvarname, expo_name, vartype, estimate, p.value) |> rename(estimate_unadjusted=estimate, p.value_unadjusted=p.value), by=c("evarname", "pvarname", "expo_name", "vartype"))
mvr2 <- tbl(con, 'mvr2') |> mutate(mv = mve_rsq-base_rsq)
pe <- tbl(con, "pe")
glanced <- tbl(con, "glanced")
variable_domain <- tbl(con, "variable_domain")

expos <- pe |> filter(term %like% 'expo%')
expos_wide <- ep_assoc_summary_across_models(expos, glanced)
expos_wide <- expos_wide |> left_join(varnames, by=c("evarname"="Variable.Name", "exposure_table_name"="Data.File.Name"))
expos_wide <- expos_wide |> left_join(varnames |> select(Variable.Name, Data.File.Name, Variable.Description, Data.File.Description), 
                                      by=c("pvarname"="Variable.Name", "phenotype_table_name"="Data.File.Name"))

expos_wide <- expos_wide |> select(-Use.Constraints) |> rename(e_data_file_desc=Data.File.Description.x, p_data_file_desc=Data.File.Description.y,
                                                               e_variable_description=Variable.Description.x, 
                                                               p_variable_description=Variable.Description.y
                                                               )
expos_wide <- expos_wide |> collect()
#expos_wide_summary <- expos_wide |> filter(term == 'expo' | term == 'expo1') |> group_by(evarname, pvarname) |> summarize(mean_adjusted_base_r2_diff = mean(rsq_adjusted_base_diff), mean_unadjusted_r2_diff=mean(rsq_adjusted_diff), total_n = sum(nobs_adjusted)) |> ungroup()

expos_wide_summary <- expos_wide |> filter(term == 'expo' | term == 'expo1' | term == 'expo2') |> group_by(evarname, pvarname) |> summarize(mean_adjusted_base_r2_diff = mean(rsq_adjusted_base_diff), mean_unadjusted_r2_diff=mean(rsq_adjusted_diff), total_n = sum(nobs_adjusted)) |> ungroup()

## `summarise()` has grouped output by 'evarname'. You can override using the
## `.groups` argument.

adjusted_meta <- adjusted_meta |> collect() |> left_join(expos_wide_summary, by=c("evarname", "pvarname"))


p_variable_domain <- variable_domain |> filter(epcf == 'p') |> collect() |> group_by(Variable.Name) |> summarise(pvardesc=first(Variable.Description),pcategory=first(category),psubcategory=first(subcategory))
e_variable_domain <- variable_domain |> filter(epcf == 'e') |> collect() |> group_by(Variable.Name) |> summarise(evardesc=first(Variable.Description),ecategory=first(category),esubcategory=first(subcategory))

adjusted_meta <- adjusted_meta |> left_join(p_variable_domain, by=c("pvarname"="Variable.Name"))
adjusted_meta <- adjusted_meta |> left_join(e_variable_domain, by=c("evarname"="Variable.Name"))

expos_wide <- expos_wide |> left_join(p_variable_domain, by=c("pvarname"="Variable.Name"))
expos_wide <- expos_wide |> left_join(e_variable_domain, by=c("evarname"="Variable.Name"))

mvr2 <- mvr2 |> collect() |> left_join(p_variable_domain, by=c("pvarname"="Variable.Name"))

Number of unique exposures and phenotypes

num_e <- length(unique(adjusted_meta$evarname))
num_p <- length(unique(adjusted_meta$pvarname))
num_e

## [1] 859

num_p

## [1] 319

num_e * num_p

## [1] 274021

Number exposures and phenotypes and associations in X number of surveys

#adjusted_meta <- adjusted_meta |> unnest(glanced) |> unnest(tidied)
n_obss <- sort(unique(adjusted_meta$nobs))

num_tests <- map_df(n_obss, function(n) {
  n_e <- adjusted_meta |> filter(nobs == n) |> pull(evarname) |> unique() |> length()
  n_p <- adjusted_meta |> filter(nobs == n) |> pull(pvarname) |> unique() |> length()
  nn <- nrow(adjusted_meta |> filter(nobs == n))
  tibble(n_expos=n_e, n_phenos=n_p, n_pxe=nn)
})

num_tests  |> mutate(n_surveys=n_obss) |> gt()

n_expos	n_phenos	n_pxe	n_surveys
853	319	58906	1
650	278	38251	2
519	241	21600	3
411	224	36138	4
319	216	9160	5
326	82	8794	6
217	75	1827	7
196	64	9423	8
38	60	1051	9
27	54	1573	10

Keep number of surveys is greater than 2

adjusted_meta_2 <- adjusted_meta |> filter(nobs >= 2)
n_evars <- length(unique(adjusted_meta_2$evarname))
n_pvars <- length(unique(adjusted_meta_2$pvarname))

Sample sizes within and across all surveys

sample_size_per_pair <- expos_wide |> filter(term == 'expo' | term== 'expo1') |> group_by(evarname, pvarname) |> summarize(total_n=sum(nobs_adjusted), n_surveys=n(), median_n=median(nobs_adjusted))

## `summarise()` has grouped output by 'evarname'. You can override using the
## `.groups` argument.

Summary of the summary stats

adjusted_meta_2 <- adjusted_meta_2 |> ungroup() |>  mutate(pval_BY=p.adjust(p.value, method="BY"), pvalue_bonferroni=p.adjust(p.value, method="bonferroni"))
adjusted_meta_2 <- adjusted_meta_2 |> mutate(sig_levels = case_when(
  pvalue_bonferroni < 0.05 ~ 'Bonf.<0.05',
  pval_BY < 0.05 ~ 'BY<0.05',
  TRUE ~ '> BY & Bonf.'
))

bonf_thresh <- 0.05/nrow(adjusted_meta_2)
quantile(adjusted_meta_2$p.value, probs=c(0.01, .05, .1, .2, .3, .4, .5, .6, .7, .8, .9, .95, .99), na.rm = T)

##           1%           5%          10%          20%          30%          40% 
## 2.553508e-27 3.398496e-10 3.332952e-06 1.761850e-03 2.141903e-02 8.095805e-02 
##          50%          60%          70%          80%          90%          95% 
## 1.826276e-01 3.179099e-01 4.734236e-01 6.422462e-01 8.190745e-01 9.093611e-01 
##          99% 
## 9.816266e-01

quantile(adjusted_meta_2$estimate, probs=c(0.01, .05, .1, .2, .3, .4, .5, .6, .7, .8, .9, .95, .99), na.rm = T)

##           1%           5%          10%          20%          30%          40% 
## -0.203919498 -0.075944258 -0.046650432 -0.023478563 -0.012016041 -0.003777513 
##          50%          60%          70%          80%          90%          95% 
##  0.003753614  0.011768619  0.021248856  0.033866834  0.057029843  0.090201321 
##          99% 
##  0.287943341

sum(adjusted_meta_2$pvalue_bonferroni < 0.05)/nrow(adjusted_meta_2) 

## [1] 0.08294671

adjusted_meta_2 |> group_by(sig_levels) |> count()

## # A tibble: 3 × 2
## # Groups:   sig_levels [3]
##   sig_levels        n
##   <chr>         <int>
## 1 > BY & Bonf. 105092
## 2 Bonf.<0.05    10602
## 3 BY<0.05       12123

adjusted_meta_2 |> filter(sig_levels == 'BY<0.05') |> arrange(-p.value) |> head()

## # A tibble: 6 × 39
##   evarname pvarname term  type  estimate std.e…¹ stati…² p.value i.squ…³ h.squ…⁴
##   <chr>    <chr>    <chr> <chr>    <dbl>   <dbl>   <dbl>   <dbl>   <dbl>   <dbl>
## 1 URXDMA   LBXPLTSI over… summ… -0.0343  0.0101    -3.38 7.20e-4     0      1   
## 2 DR2TZINC DXXRAFAT over… summ… -0.0340  0.0101    -3.38 7.20e-4     0      1   
## 3 DR1TP182 DXXRLA   over… summ…  0.0292  0.00863    3.38 7.20e-4     0      1   
## 4 LBXV1D   LBXSATSI over… summ… -0.00719 0.00213   -3.38 7.20e-4    91.5   11.8 
## 5 DSQTFA   DXDSTTOT over… summ…  0.0483  0.0143     3.38 7.20e-4     0      1   
## 6 DR1TNIAC LBXSIR   over… summ…  0.0273  0.00806    3.38 7.20e-4    19.0    1.23
## # … with 29 more variables: tau.squared <dbl>, tau.squared.se <dbl>,
## #   cochran.qe <dbl>, p.value.cochran.qe <dbl>, cochran.qm <dbl>,
## #   p.value.cochran.qm <dbl>, df.residual <int>, logLik <dbl>, deviance <dbl>,
## #   AIC <dbl>, BIC <dbl>, AICc <dbl>, nobs <int>, vartype <chr>,
## #   expo_name <chr>, estimate_unadjusted <dbl>, p.value_unadjusted <dbl>,
## #   mean_adjusted_base_r2_diff <dbl>, mean_unadjusted_r2_diff <dbl>,
## #   total_n <int>, pvardesc <chr>, pcategory <chr>, psubcategory <chr>, …
## # ℹ Use `colnames()` to see all variable names

adjusted_meta_2 |> group_by(sig_levels) |> summarize(r2_25=quantile(mean_adjusted_base_r2_diff, probs=.25, na.rm = T),
                                                     r2_50=quantile(mean_adjusted_base_r2_diff, probs=.5, na.rm = T),
                                                     r2_75=quantile(mean_adjusted_base_r2_diff, probs=.75, na.rm = T),
                                                      r2_100=max(mean_adjusted_base_r2_diff, na.rm = T))

## # A tibble: 3 × 5
##   sig_levels      r2_25    r2_50   r2_75 r2_100
##   <chr>           <dbl>    <dbl>   <dbl>  <dbl>
## 1 > BY & Bonf. 0.000357 0.000763 0.00162 0.112 
## 2 Bonf.<0.05   0.00225  0.00425  0.00842 0.547 
## 3 BY<0.05      0.00114  0.00196  0.00371 0.0785

adjusted_meta_2 |> filter(sig_levels == 'Bonf.<0.05') |> arrange(-mean_adjusted_base_r2_diff)

## # A tibble: 10,602 × 39
##    evarn…¹ pvarn…² term  type  estim…³ std.e…⁴ stati…⁵   p.value i.squ…⁶ h.squ…⁷
##    <chr>   <chr>   <chr> <chr>   <dbl>   <dbl>   <dbl>     <dbl>   <dbl>   <dbl>
##  1 URXDHB  URXUCR  over… summ…   0.766  0.0191    40.1 0             0      1   
##  2 URXPHG  URXUCR  over… summ…   0.673  0.0214    31.4 4.42e-216     0      1   
##  3 URXUCS  URXUCR  over… summ…   0.647  0.0214    30.2 1.39e-200     0      1   
##  4 URXMAD  URXUCR  over… summ…   0.662  0.0235    28.1 4.19e-174     0      1   
##  5 URXMBP  URXUCR  over… summ…   0.642  0.0180    35.7 1.97e-279     0      1   
##  6 URXUMO  URXUCR  over… summ…   0.633  0.0184    34.3 3.01e-258     0      1   
##  7 LBXOTT  LBXSTR  over… summ…   0.633  0.0391    16.2 7.28e- 59    39.9    1.66
##  8 LBXOTT  LBXTR   over… summ…   0.631  0.0392    16.1 3.09e- 58    38.2    1.62
##  9 URXMIB  URXUCR  over… summ…   0.602  0.0202    29.9 3.58e-196     0      1   
## 10 URXUTL  URXUCR  over… summ…   0.585  0.0158    36.9 8.25e-299     0      1   
## # … with 10,592 more rows, 29 more variables: tau.squared <dbl>,
## #   tau.squared.se <dbl>, cochran.qe <dbl>, p.value.cochran.qe <dbl>,
## #   cochran.qm <dbl>, p.value.cochran.qm <dbl>, df.residual <int>,
## #   logLik <dbl>, deviance <dbl>, AIC <dbl>, BIC <dbl>, AICc <dbl>, nobs <int>,
## #   vartype <chr>, expo_name <chr>, estimate_unadjusted <dbl>,
## #   p.value_unadjusted <dbl>, mean_adjusted_base_r2_diff <dbl>,
## #   mean_unadjusted_r2_diff <dbl>, total_n <int>, pvardesc <chr>, …
## # ℹ Use `print(n = ...)` to see more rows, and `colnames()` to see all variable names

## qqplot
pval_qq <- data.frame(observed = sort(adjusted_meta_2$p.value), expected = (1:nrow(adjusted_meta_2))/nrow(adjusted_meta_2))
qq_p <- ggplot(pval_qq, aes(-log10(expected), -log10(observed)))
qq_p <- qq_p + geom_point()
##

p <- ggplot(pval_qq, aes(observed))
p <- p + geom_histogram(bins=100) + theme_bw()
p <- p + geom_hline(yintercept = 1, color='blue')
p <- p + xlab("P-E association pvalue")
p_hist <- p
p_hist

p <- ggplot(pval_qq |> filter(observed < 1e-3), aes(-log10(observed)))
p <- p + geom_histogram(bins=200) + theme_bw() + scale_x_continuous(limits=c(0, 100))
p <- p + geom_hline(yintercept = 1, color='blue')
p <- p + xlab("P-E association pvalue")
p_hist <- p
p_hist

## Warning: Removed 122 rows containing non-finite values (stat_bin).

## Warning: Removed 2 rows containing missing values (geom_bar).

p <- ggplot(adjusted_meta_2, aes(p.value))
p <- p + geom_density() + theme_bw() + facet_grid(~ecategory)
p1 <- p + xlab("P-E association pvalue")


p <- ggplot(adjusted_meta_2, aes(p.value))
p <- p + geom_density() + theme_bw() + facet_grid(~pcategory)
p2 <- p + xlab("P-E association pvalue")

plot_grid(p1, p2, nrow=2, labels=c("A", "B"))

zoom in the distribution

p_plot <- adjusted_meta_2 |> select(ecategory, pcategory, p.value) 

p <- ggplot(p_plot,aes(p.value))
p <- p + geom_histogram(aes(y=..density..)) + geom_density()  + theme_bw() + facet_grid(~ecategory)
p1 <- p + xlab("P-E association pvalue")


p <- ggplot(p_plot, aes(p.value))
p <- p + geom_histogram(aes(y=..density..)) + geom_density() + theme_bw() + facet_grid(~pcategory)
p2 <- p + xlab("P-E association pvalue")

plot_grid(p1, p2, nrow=2, labels=c("A", "B"))

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Association Size vs. -log10(pvalue)

e_summary <- adjusted_meta_2 |> group_by(evarname) |> arrange(pvalue_bonferroni) |>  
  summarize(mean_r2=mean(mean_adjusted_base_r2_diff, na.rm=T),  mean_estimate=mean(abs(estimate), na.rm=T), 
            median_pvalue=median(p.value, na.rm=T), n_sig=sum(pvalue_bonferroni < 0.05, na.rm=T), 
            n_tests=sum(!is.na(pvalue_bonferroni)),  median_i.squared=median(i.squared, na.rm=T),
            max_r2=first(mean_adjusted_base_r2_diff), max_pvarname=first(pvarname) , max_estimate=first(estimate), max_p.value=first(p.value)) |> mutate(n_sig_pct=n_sig/n_tests)


p_summary <- adjusted_meta_2 |> group_by(pvarname) |> arrange(pvalue_bonferroni) |> 
  summarize(mean_r2=mean(mean_adjusted_base_r2_diff, na.rm=T), mean_estimate=mean(abs(estimate), na.rm=T),
            median_pvalue=median(p.value, na.rm=T), n_sig=sum(pvalue_bonferroni < 0.05, na.rm=T), 
            n_tests=sum(!is.na(pvalue_bonferroni)),  median_i.squared=median(i.squared, na.rm=T),
            max_r2=first(mean_adjusted_base_r2_diff), max_evarname=first(evarname) , max_estimate=first(estimate), max_p.value=first(p.value)) |> mutate(n_sig_pct=n_sig/n_tests)

## deeper summary by group

p_group_summary <- adjusted_meta_2 |> unite(p_scategory, c(pcategory, psubcategory)) |> group_by(p_scategory) |> arrange(pvalue_bonferroni) |>  
  summarize(mean_r2=mean(mean_adjusted_base_r2_diff, na.rm=T),  mean_estimate=mean(abs(estimate), na.rm=T), 
            median_pvalue=median(p.value, na.rm=T), n_sig=sum(pvalue_bonferroni < 0.05, na.rm=T), 
            n_tests=sum(!is.na(pvalue_bonferroni)),  median_i.squared=median(i.squared, na.rm=T),
            max_r2=first(mean_adjusted_base_r2_diff), max_evarname=first(evarname) , max_estimate=first(estimate), max_p.value=first(p.value)) |> mutate(n_sig_pct=n_sig/n_tests)


e_group_summary <- adjusted_meta_2 |> unite(e_scategory, c(ecategory, esubcategory)) |> group_by(e_scategory) |> arrange(pvalue_bonferroni) |>  
  summarize(mean_r2=mean(mean_adjusted_base_r2_diff, na.rm=T),  
            mean_abs_estimate=mean(abs(estimate), na.rm=T),
            mean_estimate=mean((estimate), na.rm=T),
            median_pvalue=median(p.value, na.rm=T), 
            n_sig=sum(pvalue_bonferroni < 0.05, na.rm=T), 
            n_tests=sum(!is.na(pvalue_bonferroni)),  
            median_i.squared=median(i.squared, na.rm=T),
            max_r2=first(mean_adjusted_base_r2_diff), 
            max_pvarname=first(pvarname), 
            max_estimate=first(estimate), max_p.value=first(p.value)) |> mutate(n_sig_pct=n_sig/n_tests)

Tables 1 and 2

## 
e_group_summary <- adjusted_meta_2 |> filter(expo_name %in% c('expo', 'expo1', 'expo2', 'expo3')) |> filter(mean_adjusted_base_r2_diff <= .1) |> unite(e_scategory, c(ecategory, esubcategory)) |> group_by(e_scategory) |> arrange(pvalue_bonferroni,mean_adjusted_base_r2_diff) |>  
  summarize(
            #median_r2=mean(mean_adjusted_base_r2_diff, na.rm=T),  
            #median_abs_estimate=median(abs(estimate), na.rm=T),
            n_sig=sum(pvalue_bonferroni < 0.05, na.rm=T), 
            n_tests=sum(!is.na(pvalue_bonferroni)),  
            #median_i.squared=median(i.squared, na.rm=T),
            max_r2=first(mean_adjusted_base_r2_diff), 
            max_termname=remove_units_from_string(first(expo_name)),
            max_pvarname=remove_units_from_string(first(pvardesc)),
            max_evarname=remove_units_from_string(first(evardesc)),
            max_estimate=first(estimate), max_p.value=first(p.value), max_i.squared=first(i.squared)) |> mutate(n_sig_pct=n_sig/n_tests)

e_bonf_group_summary <- adjusted_meta_2 |> filter(sig_levels == 'Bonf.<0.05') |> 
  unite(e_scategory, c(ecategory, esubcategory)) |>
  mutate(sgn=ifelse(sign(estimate) <0, 'neg', 'pos')) |> group_by(e_scategory, sgn) |> arrange(pvalue_bonferroni) |>  
  summarize(median_bonf_r2=median(mean_adjusted_base_r2_diff, na.rm=T),  
            q25_bonf_estimate=quantile(estimate, probs=.25, na.rm=T),
            median_bonf_estimate=median((estimate), na.rm=T),
            q75_bonf_estimate=quantile(estimate, probs=.75, na.rm=T)
  )

## `summarise()` has grouped output by 'e_scategory'. You can override using the
## `.groups` argument.

e_bonf_group_summary <- e_bonf_group_summary |> pivot_wider(names_from=sgn, values_from=(c(median_bonf_r2, q25_bonf_estimate, median_bonf_estimate, q75_bonf_estimate)))

## merge
e_group_summary <- e_group_summary |> left_join(e_bonf_group_summary, by='e_scategory')
e_group_summary <- e_group_summary |> separate(col=e_scategory, into=c("escategory", "esubcategory"), sep="_")
e_group_summary <- e_group_summary |> filter(n_sig > 1)

e_group_summary <- e_group_summary |> select(escategory, esubcategory, n_tests, n_sig_pct, median_bonf_r2_neg, median_bonf_r2_pos, q25_bonf_estimate_neg, median_bonf_estimate_neg, q75_bonf_estimate_neg, q25_bonf_estimate_pos, median_bonf_estimate_pos, q75_bonf_estimate_pos, 
max_pvarname, max_evarname, max_estimate, max_r2, max_p.value, max_i.squared)

e_group_summary<- e_group_summary |> gt() |> 
  fmt_number(columns = c(starts_with("q"), starts_with("med")), decimals = 3) |>
  fmt_number(columns = c(max_r2, max_estimate), decimals = 3) |>
  fmt_number(columns = c(n_sig_pct), decimals = 2) |>
  fmt_number(columns = c(max_i.squared), decimals = 0) |>
  fmt_scientific(columns = max_p.value, decimals = 0)

e_group_summary <- e_group_summary |> 
  tab_spanner(label = "Effect Sizes < 0",columns = ends_with("neg")) |>
  tab_spanner(label = "Effect Sizes > 0",columns = ends_with("pos"))

e_group_summary <- e_group_summary |> 
  tab_spanner(label = "Example P-E Association (lowest p.value for E domain)",columns = starts_with("max"))

e_group_summary |> gtsave('./e_group_summary.html')
e_group_summary

escategory	esubcategory	n_tests	n_sig_pct	Effect Sizes < 0				Effect Sizes > 0				Example P-E Association (lowest p.value for E domain)
				median_bonf_r2_neg	q25_bonf_estimate_neg	median_bonf_estimate_neg	q75_bonf_estimate_neg	median_bonf_r2_pos	q25_bonf_estimate_pos	median_bonf_estimate_pos	q75_bonf_estimate_pos	max_pvarname	max_evarname	max_estimate	max_r2	max_p.value	max_i.squared
allergy	NA	464	0.01	NA	NA	NA	NA	0.001	0.321	0.331	0.650	Ward's triangle bone mineral content	Tissue transglutaminase	0.650	0.002	1 × 10−10	0
housing	NA	556	0.01	0.001	−0.603	−0.472	−0.258	0.001	0.174	0.275	0.384	Alkaline phosphotase	Is this {mobile home/house/apartment} owned, being bought, rented, or occup	−0.103	0.002	3 × 10−17	27
income	NA	513	0.01	0.005	−0.464	−0.421	−0.390	0.007	0.203	0.211	0.218	Left Arm Total	Do {you/NAMES OF OTHER FAMILY/you and NAMES OF FAMILY MEMBERS} have more th	0.194	0.006	5 × 10−8	0
infection	NA	7015	0.04	0.002	−1.007	−0.450	−0.196	0.002	0.115	0.239	0.514	Subtotal Lean excl BMC	Hepatitis D	0.660	0.000	2 × 10−86	0
nutrients	dietary biomarker	8301	0.16	0.010	−0.153	−0.104	−0.076	0.009	0.078	0.105	0.169	Total cholesterol	g-tocopherol	0.263	0.064	1 × 10−233	0
nutrients	dietary interview	37852	0.12	0.004	−0.076	−0.061	−0.049	0.002	0.041	0.048	0.060	Direct HDL-Cholesterol	Alcohol	0.201	0.037	3 × 10−145	0
nutrients	supplements	14619	0.01	0.004	−0.120	−0.076	−0.057	0.012	0.086	0.117	0.144	RBC folate	Any Dietary Supplements taken in the past 24 hour?	−0.511	0.050	2 × 10−178	0
physical activity	NA	4126	0.18	0.005	−0.105	−0.068	−0.054	0.004	0.049	0.064	0.096	Left Leg Percent Fat	Vigorous Leisure Per Week	−0.155	0.020	5 × 10−91	2
pollutant	amide	456	0.22	0.012	−0.138	−0.112	−0.081	0.005	0.064	0.078	0.108	White blood cell count	Acrylamide	0.142	0.018	4 × 10−45	0
pollutant	amine	1233	0.05	0.007	−0.106	−0.086	−0.082	0.013	0.098	0.116	0.162	Creatinine, urine	4-Aminobiphenyl, urine	0.322	0.095	3 × 10−104	0
pollutant	deet	712	0.00	0.003	−0.054	−0.054	−0.054	0.006	0.074	0.074	0.074	Uric acid	DEET acid	0.074	0.006	4 × 10−8	0
pollutant	diakyl	2285	0.05	0.002	−0.066	−0.042	−0.021	0.004	0.053	0.063	0.077	Blood urea nitrogen	Urinary Perchlorate	0.126	0.016	4 × 10−74	1
pollutant	heavy metals	7450	0.11	0.006	−0.108	−0.086	−0.067	0.006	0.066	0.081	0.097	Waist Circumference	Lead	−0.164	0.018	2 × 10−174	0
pollutant	hydrocarbon	2388	0.11	0.008	−0.116	−0.091	−0.074	0.007	0.071	0.089	0.105	White blood cell count	2-fluorene	0.172	0.028	1 × 10−57	6
pollutant	organochlorine	5231	0.04	0.030	−0.345	−0.264	−0.151	0.022	0.131	0.181	0.271	Arm Circumference	PCB194	−0.374	0.050	1 × 10−45	0
pollutant	organophosphate	472	0.04	0.000	−0.034	−0.018	−0.014	0.000	0.011	0.011	0.012	Segmented neutrophils percent	Dimethoate	−0.028	0.000	3 × 10−94	0
pollutant	phenols	2887	0.03	0.004	−0.073	−0.065	−0.057	0.006	0.066	0.082	0.093	Body Mass Index	Ethyl paraben	−0.110	0.012	5 × 10−36	0
pollutant	phthalates	4390	0.11	0.007	−0.095	−0.081	−0.068	0.008	0.068	0.085	0.100	Osmolality	Mono-	0.142	0.019	4 × 10−54	0
pollutant	polyfluoro	3111	0.06	0.010	−0.130	−0.106	−0.078	0.012	0.091	0.115	0.158	Albumin	Perfluorohexane sulfonic acid	0.157	0.023	9 × 10−40	6
pollutant	priority pesticide	2574	0.03	0.000	−0.024	−0.014	−0.011	0.001	0.013	0.034	0.078	Arm Circumference	2,4,5-T (ug/L) result	−0.026	0.001	2 × 10−249	0
pollutant	VOC	12975	0.06	0.006	−0.095	−0.079	−0.050	0.007	0.056	0.080	0.109	Baseline 1st Test Spirometry, Forced Vital Capacity, in mL.	Blood Carbon Tetrachloride Result	−0.015	0.000	1 × 10−96	0
smoking	smoking behavior	2717	0.04	0.007	−0.254	−0.205	−0.150	0.008	0.106	0.162	0.277	Mean cell volume	Smoking: current, ever, never	0.474	0.031	2 × 10−194	21
smoking	smoking biomarker	1393	0.11	0.004	−0.088	−0.072	−0.057	0.010	0.065	0.108	0.162	Segmented neutrophils num	Cotinine	0.165	0.023	6 × 10−211	0

## 
p_group_summary <- adjusted_meta_2 |> filter(expo_name %in% c('expo', 'expo1', 'expo2', 'expo3')) |> filter(mean_adjusted_base_r2_diff <= .1) |> unite(p_scategory, c(pcategory, psubcategory)) |> group_by(p_scategory) |> arrange(pvalue_bonferroni,mean_adjusted_base_r2_diff) |>  
  summarize(
            #median_r2=mean(mean_adjusted_base_r2_diff, na.rm=T),  
            #median_abs_estimate=median(abs(estimate), na.rm=T),
            n_sig=sum(pvalue_bonferroni < 0.05, na.rm=T), 
            n_tests=sum(!is.na(pvalue_bonferroni)),  
            #median_i.squared=median(i.squared, na.rm=T),
            max_r2=first(mean_adjusted_base_r2_diff), 
            max_termname=remove_units_from_string(first(expo_name)),
            max_pvarname=remove_units_from_string(first(pvardesc)),
            max_evarname=remove_units_from_string(first(evardesc)),
            max_estimate=first(estimate), max_p.value=first(p.value), max_i.squared=first(i.squared)) |> mutate(n_sig_pct=n_sig/n_tests)

p_bonf_group_summary <- adjusted_meta_2 |> filter(sig_levels == 'Bonf.<0.05') |> 
  unite(p_scategory, c(pcategory, psubcategory)) |>
  mutate(sgn=ifelse(sign(estimate) <0, 'neg', 'pos')) |> group_by(p_scategory, sgn) |> arrange(pvalue_bonferroni) |>  
  summarize(median_bonf_r2=median(mean_adjusted_base_r2_diff, na.rm=T),  
            q25_bonf_estimate=quantile(estimate, probs=.25, na.rm=T),
            median_bonf_estimate=median((estimate), na.rm=T),
            q75_bonf_estimate=quantile(estimate, probs=.75, na.rm=T)
  )

## `summarise()` has grouped output by 'p_scategory'. You can override using the
## `.groups` argument.

p_bonf_group_summary <- p_bonf_group_summary |> pivot_wider(names_from=sgn, values_from=(c(median_bonf_r2, q25_bonf_estimate, median_bonf_estimate, q75_bonf_estimate)))

## merge
p_group_summary <- p_group_summary |> left_join(p_bonf_group_summary, by='p_scategory')
p_group_summary <- p_group_summary |> separate(col=p_scategory, into=c("pscategory", "psubcategory"), sep="_")
p_group_summary <- p_group_summary |> filter(n_sig > 1)

p_group_summary <- p_group_summary |> select(pscategory, psubcategory, n_tests, n_sig_pct, q25_bonf_estimate_neg, median_bonf_estimate_neg, q75_bonf_estimate_neg, q25_bonf_estimate_pos, median_bonf_estimate_pos, q75_bonf_estimate_pos, 
max_pvarname, max_evarname, max_estimate, max_r2, max_p.value, max_i.squared)

p_group_summary<- p_group_summary |> gt() |> 
  fmt_number(columns = c(starts_with("q"), starts_with("med")), decimals = 3) |>
  fmt_number(columns = c(max_r2, max_estimate), decimals = 3) |>
  fmt_number(columns = c(n_sig_pct), decimals = 2) |>
  fmt_number(columns = c(max_i.squared), decimals = 0) |>
  fmt_scientific(columns = max_p.value, decimals = 0)

p_group_summary <- p_group_summary |> 
  tab_spanner(label = "Effect Sizes < 0",columns = ends_with("neg")) |>
  tab_spanner(label = "Effect Sizes > 0",columns = ends_with("pos"))

p_group_summary <- p_group_summary |> 
  tab_spanner(label = "Example P-E Association (lowest p.value for E domain)",columns = starts_with("max"))

p_group_summary |> gtsave('./p_group_summary.html')
p_group_summary

pscategory	psubcategory	n_tests	n_sig_pct	Effect Sizes < 0			Effect Sizes > 0			Example P-E Association (lowest p.value for E domain)
				q25_bonf_estimate_neg	median_bonf_estimate_neg	q75_bonf_estimate_neg	q25_bonf_estimate_pos	median_bonf_estimate_pos	q75_bonf_estimate_pos	max_pvarname	max_evarname	max_estimate	max_r2	max_p.value	max_i.squared
aging	NA	362	0.01	NA	NA	NA	0.183	0.310	0.437	Mean T/S ratio	Platinum, urine	0.056	0.002	7 × 10−20	6
anthropometric	dexa	57728	0.09	−0.100	−0.076	−0.061	0.044	0.056	0.082	Android fat mass	25-hydroxyvitamin D3	−0.217	0.036	3 × 10−133	0
anthropometric	NA	6650	0.18	−0.121	−0.075	−0.050	0.041	0.060	0.088	Arm Circumference	2,4,5-T (ug/L) result	−0.026	0.001	2 × 10−249	0
biochemistry	bone	1381	0.08	−0.081	−0.072	−0.041	0.045	0.066	0.099	Total calcium	Retinol	0.230	0.044	6 × 10−185	0
biochemistry	electrolyte	3308	0.09	−0.093	−0.073	−0.059	0.053	0.089	0.112	Chloride	Molybdenum, urine	0.176	0.029	6 × 10−84	0
biochemistry	hormone	4759	0.02	−0.107	−0.080	−0.071	0.049	0.078	0.105	Parathyroid Hormone(Elecys method) pg/mL	Vitamin D	−0.269	0.058	2 × 10−71	14
biochemistry	immunity	662	0.10	−0.078	−0.059	−0.042	0.064	0.075	0.098	Globulin	Metsulfuron methyl	−0.028	0.001	6 × 10−51	0
biochemistry	inflammation	708	0.09	−0.084	−0.064	−0.054	0.060	0.082	0.111	C-reactive protein	cis-b-carotene	−0.159	0.022	1 × 10−66	0
biochemistry	injury	591	0.03	−0.074	−0.066	−0.051	0.050	0.060	0.080	Lactate dehydrogenase LDH	Smoking: current, ever, never	−0.169	0.004	8 × 10−42	0
biochemistry	kidney	2897	0.17	−0.085	−0.060	−0.046	0.068	0.104	0.324	Uric acid	Retinol	0.207	0.034	1 × 10−157	19
biochemistry	lipids	4039	0.06	−0.099	−0.069	−0.055	0.077	0.199	0.311	Total cholesterol	g-tocopherol	0.263	0.064	1 × 10−233	0
biochemistry	liver	3170	0.06	−0.089	−0.063	−0.042	0.052	0.066	0.088	Alkaline phosphotase	Pyridoxal 5'-phosphate	−0.115	0.013	8 × 10−84	0
biochemistry	liver/kidney	1994	0.10	−0.089	−0.068	−0.046	0.049	0.073	0.113	Albumin	Pyridoxal 5'-phosphate	0.252	0.060	1 × 10−84	52
biochemistry	metabolic	3140	0.07	−0.106	−0.073	−0.052	0.072	0.127	0.174	Glycohemoglobin	Nicosulfuron	0.016	0.000	2 × 10−151	0
biochemistry	nutritional status	5042	0.11	−0.116	−0.091	−0.074	0.074	0.092	0.122	Folate, RBC	Vitamin B12	0.245	0.057	7 × 10−200	0
biochemistry	psa	1179	0.01	−0.088	−0.086	−0.072	0.153	0.153	0.153	Prostate specific antigen ratio	Cadmium	−0.086	0.007	3 × 10−9	0
blood pressure	NA	2952	0.05	−0.096	−0.065	−0.045	0.061	0.079	0.104	60 sec. pulse (30 sec. pulse * 2):	trans-b-carotene	−0.144	0.018	1 × 10−67	0
blood	iron	3549	0.04	−0.152	−0.109	−0.075	0.062	0.091	0.143	Iron, refigerated	g-tocopherol	−0.151	0.021	7 × 10−110	0
blood	NA	12763	0.06	−0.080	−0.051	−0.039	0.064	0.095	0.129	Segmented neutrophils num	Cotinine	0.165	0.023	6 × 10−211	0
fitness	NA	666	0.06	−0.124	−0.088	−0.069	0.090	0.113	0.129	Predicted maximal oxygen uptake	g-tocopherol	−0.156	0.023	5 × 10−76	0
lung	exhaled NO	390	0.04	−0.196	−0.171	−0.140	0.211	0.373	0.536	Mean of two reproducible FENO measurements, in parts per billion	Smoking: current, ever, never	−0.524	0.040	2 × 10−86	0
lung	lung function	782	0.09	−0.143	−0.069	−0.046	0.041	0.047	0.059	Baseline 1st Test Spirometry, Forced Vital Capacity, in mL.	Blood Carbon Tetrachloride Result	−0.015	0.000	1 × 10−96	0
microbiome	NA	5520	0.01	−0.076	−0.070	−0.034	0.026	0.104	0.117	RSV_ObservedOTUsSD	Herpes Simplex Virus II	−1.065	0.001	4 × 10−29	0

adjusted_meta_2 <- adjusted_meta_2 |> mutate(p_cap = ifelse(p.value < 1e-30, 1e-30, p.value))
p <- ggplot(adjusted_meta_2 |> filter(ecategory != 'autoantibody'), aes(estimate, -log10(p_cap)))
p <- p + geom_point(shape='.') + scale_x_continuous(limits=c(-1, 1))
p <- p + facet_grid(ecategory ~ .) + scale_color_npg()
p <- p + geom_hline(yintercept = -log10(.05/nrow(adjusted_meta_2)), color='lightblue')
p <- p + theme_minimal() + theme(legend.position = "none") +ylab('p.value') + xlab("estimate")
p1 <- p

## uBiome only
p <- ggplot(adjusted_meta_2 |> filter(pcategory == 'microbiome'), aes(estimate, -log10(p_cap)))
p <- p + geom_point(shape='.') + scale_x_continuous(limits=c(-1, 1))
p <- p + facet_grid(ecategory ~ .) + scale_color_npg()
p <- p + geom_hline(yintercept = -log10(.05/nrow(adjusted_meta_2)), color='lightblue')
p <- p + theme_minimal() + theme(legend.position = "none") +ylab('p.value') + xlab("estimate")


p <- ggplot(adjusted_meta_2, aes(estimate, -log10(p_cap)))
p <- p + geom_point(shape='.') + scale_x_continuous(limits=c(-1, 1))
p <- p + facet_grid(pcategory ~ .) + scale_color_npg()
p <- p + geom_hline(yintercept = -log10(.05/nrow(adjusted_meta_2)), color='lightblue')
p <- p + theme_minimal() + theme(legend.position = "none") +ylab('p.value') + xlab("estimate")
p2 <-p

plot_grid(p1, p2, ncol=2, labels=c("A", "B"))

## Warning: Removed 73 rows containing missing values (geom_point).
## Removed 73 rows containing missing values (geom_point).

# library(circlize)
# # create circos plot visualizing 5 phenotype associations from adjusted_meta_2
# pvarnames_to_circos <- c("BMXBMI", "MSYSTOLIC","LBXGLU", "TeloMean", "LBXTC")
# pcircos <- adjusted_meta_2 |> filter(pvarname %in% pvarnames_to_circos) |> filter(vartype == 'continuous') |> select(c("evarname", "ecategory", "esubcategory", "evardesc", "pvarname", "pval_BY", "estimate"))
# 
# pcircos <- pcircos |> unite(ecategory_sub, ecategory, esubcategory, sep="_")
# 
# evarname_per_sub <- pcircos |> group_by(ecategory_sub, evarname) |> count() |> ungroup() |> select(-n)
# #m_per_ecategory <- uevarname_per_sub |> ungroup() |> select(-n) |> group_by(ecategory_sub) |> count()
# 
# evarname_per_sub <- evarname_per_sub |> group_by(ecategory_sub) |> mutate(index = match(evarname, unique(evarname)))
# #m_per_ecategory <- m_per_ecategory |> filter(ecategory_sub != 'pollutant_pyrethoid',
# #                                             ecategory_sub != 'allergy_NA',
# #                                             ecategory_sub != 'pollutant_amide',
# #                                             )
# pcircos_2 <- pcircos |> left_join(evarname_per_sub |> ungroup()) # this will give an index to each
# m_per_ecategory <- evarname_per_sub |> count()
# 
# m_per_ecategory <- m_per_ecategory |> filter(n >= 10)
# 
# xlims <- matrix(nrow=length(m_per_ecategory$ecategory_sub), data=0, ncol=1) |> cbind(matrix(nrow=length(m_per_ecategory$ecategory_sub), data=m_per_ecategory$n, ncol=1))
# 
# #circos.par(cell.padding = c(0.02, 0, 0.02, 0))
# pcircos_3 <- pcircos_2 |> filter(ecategory_sub %in% m_per_ecategory$ecategory_sub)
# pcircos_3 <- pcircos_3 |> select(evarname, pvarname, estimate) |>  complete(evarname, pvarname)
# estimate_wide_tbl <- pcircos_3 |> pivot_wider(names_from = pvarname, values_from = estimate) |> left_join(pcircos_2 |> select(evarname, ecategory_sub) |> unique())
# estimate_wide_matr <- estimate_wide_tbl |> select(-evarname, -ecategory_sub) |> as.matrix()
# 
# exposome_domains <- factor(estimate_wide_tbl$ecategory_sub)
# col_fun1 = colorRamp2(c(-.1, 0, .1), c("blue", "white", "red"))
# circos.clear()
# circos.par(gap.degree = 5)
# circos.heatmap.initialize(estimate_wide_matr,split=exposome_domains)
# circos.heatmap(estimate_wide_matr[,1], col=col_fun1)
# circos.heatmap(estimate_wide_matr[,2], col=col_fun1)
# circos.heatmap(estimate_wide_matr[,3], col=col_fun1)
# circos.heatmap(estimate_wide_matr[,4], col=col_fun1)
# circos.heatmap(estimate_wide_matr[,5], col=col_fun1)

e_effect_sizes_per <- adjusted_meta_2 |> filter(sig_levels == 'Bonf.<0.05') |> group_by(ecategory, esubcategory, sign(estimate)) |> summarize(median_pvalue=median(p.value), number_signficant=n(), mean_estimate=mean((estimate))) |> arrange(-mean_estimate)

## `summarise()` has grouped output by 'ecategory', 'esubcategory'. You can
## override using the `.groups` argument.

e_effect_sizes_per <- e_effect_sizes_per |> mutate(esubcategory = ifelse(is.na(esubcategory), ecategory, esubcategory))
p <- ggplot(e_effect_sizes_per, aes(mean_estimate, -log10(median_pvalue), label=esubcategory))
p <- p + geom_point(aes(size=number_signficant)) + geom_text_repel() + geom_vline(xintercept = 0)
p <- p + theme_bw() + xlab("Average(Estimate) within exposome groups") + ylab("Median log10(pvalue)")
p <- p + theme(legend.position = "bottom")
p

## Warning: ggrepel: 14 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

p_effect_sizes_per <- adjusted_meta_2 |> filter(sig_levels == 'Bonf.<0.05') |> group_by(pcategory, psubcategory) |> summarize(mean_r2 = mean(mean_adjusted_base_r2_diff, na.rm=T))

## `summarise()` has grouped output by 'pcategory'. You can override using the
## `.groups` argument.

p <- ggplot(adjusted_meta_2 |> filter(vartype =='categorical'), aes(abs(estimate), color=sig_levels))
p <- p + stat_ecdf() + scale_x_continuous(limits=c(0, .25))
p <- p + xlab("abs(estimate)") + ylab("percentile")  + theme(legend.position="bottom") + scale_color_npg()
p

## Warning: Removed 2210 rows containing non-finite values (stat_ecdf).

p <- ggplot(adjusted_meta_2 |> filter(vartype =='continuous'), aes(abs(estimate), color=sig_levels))
p <- p + stat_ecdf() + scale_x_continuous(limits=c(0, .25))
p <- p + xlab("abs(estimate)") + ylab("percentile") + theme(legend.position="bottom") + scale_color_npg()
p

## Warning: Removed 263 rows containing non-finite values (stat_ecdf).

p <- ggplot(adjusted_meta_2, aes(abs(estimate), color=sig_levels))
p <- p + stat_ecdf() + scale_x_continuous(limits=c(0, .25))
p <- p + xlab("abs(estimate)") + ylab("percentile") + theme(legend.position="bottom") + scale_color_npg()
p

## Warning: Removed 2479 rows containing non-finite values (stat_ecdf).

ecdf_for_sig <- adjusted_meta_2 |> filter(sig_levels == 'Bonf.<0.05') |> pull(mean_adjusted_base_r2_diff) |> ecdf()
ecdf_for_non_sig <- adjusted_meta_2 |> filter(sig_levels == '> BY & Bonf.') |> pull(mean_adjusted_base_r2_diff) |> ecdf()


p_effect_sizes_per <- p_effect_sizes_per |> mutate(q = ecdf_for_sig(mean_r2), sig_levels ='Bonf.<0.05')
p_effect_sizes_per <- p_effect_sizes_per |> mutate(p_cat = ifelse(is.na(psubcategory), pcategory, psubcategory))
p <- ggplot(adjusted_meta_2, aes(mean_adjusted_base_r2_diff, color=sig_levels))
p <- p + stat_ecdf() + scale_x_continuous(limits=c(0, .05)) +scale_color_aaas()
p <- p + geom_point(data=p_effect_sizes_per, aes(x=mean_r2, y = q, color=sig_levels)) 
p <- p + geom_text_repel(data=p_effect_sizes_per, aes(x=mean_r2, y = q, color=sig_levels, label=p_cat)) 
p <- p + xlab("R^2 (adjusted-base model)") + ylab("percentile") 
p <- p + theme_bw() + theme(legend.position="bottom") 
p

## Warning: Removed 1395 rows containing non-finite values (stat_ecdf).

## Warning: ggrepel: 4 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

Association sizes for all exposures

contextualizing all exposures

adjusted_meta_2 <- adjusted_meta_2 |> mutate(evarname = fct_reorder(evarname, abs(estimate), mean))
#p <- ggplot(adjusted_meta_2, aes(x=evarname, y=estimate, group=evarname))
#p <- p + geom_density_ridges()
#p <- p + geom_boxplot() + facet_grid(~sig_levels)
#p

Replicability and Consistency

I2 across number of surveys for P-E pair

## reorder sig_levels
adjusted_meta_2 <- adjusted_meta_2 |> mutate(sig_levels=fct_relevel(sig_levels, c("> BY & Bonf.","BY<0.05", "Bonf.<0.05")))

p <- ggplot(adjusted_meta_2, aes(factor(nobs), i.squared,color=sig_levels))
p <- p + geom_boxplot() + xlab("Number of Surveys for PE Association") + theme(legend.position="bottom") + scale_color_aaas()
p <- p + theme_bw() + theme(legend.position="bottom") + ylab("i-squared")
p

i2_medians <- adjusted_meta_2 |> group_by(sig_levels) |> summarize(i2_median=median(i.squared))

Reverse meta-analysis: replicability of P-E

within_survey_pvalue_threshold <- 0.05 # should this be changed for each survey?
p_val_pe_pair <- expos_wide |> group_by(evarname, pvarname, term) |> summarize(n_pvalue_lt=sum(p.value_adjusted<within_survey_pvalue_threshold), total_n=n(), pct_pvalue_lt=n_pvalue_lt/total_n)

## `summarise()` has grouped output by 'evarname', 'pvarname'. You can override
## using the `.groups` argument.

adjusted_meta_3 <- adjusted_meta_2 |> left_join(p_val_pe_pair, by=c("evarname"="evarname", "pvarname"="pvarname", "expo_name"="term"))

Showcasing associations:

• Low pvalue, Higher R2, low I2
• Low pvalue, Higher R2, and higher I2

adjusted_meta_3 |> filter(pvalue_bonferroni < 0.05) |> nrow() 

## [1] 10602

adjusted_meta_3 |> filter(pvalue_bonferroni < 0.05, n_pvalue_lt >= 2) |> nrow()

## [1] 5048

non_het_pairs <- adjusted_meta_3 |> filter(pvalue_bonferroni < 0.05, n_pvalue_lt >= 2, i.squared < 50, mean_adjusted_base_r2_diff > .025)
het_pairs <- adjusted_meta_3 |> filter(pvalue_bonferroni < 0.05, n_pvalue_lt >= 2, i.squared > 50, mean_adjusted_base_r2_diff > .025, nobs >= 4)
#het_pairs_2 <- adjusted_meta_3 |> filter(sig_levels == 'BY<0.05', i.squared > 90, nobs >= 4) 

adjusted_meta_3 |> filter(pvalue_bonferroni < 0.05) |> group_by(ecategory) |> count()

## # A tibble: 8 × 2
## # Groups:   ecategory [8]
##   ecategory             n
##   <chr>             <int>
## 1 allergy               5
## 2 housing              40
## 3 income               10
## 4 infection           286
## 5 nutrients          6011
## 6 physical activity   740
## 7 pollutant          3236
## 8 smoking             274

## non-heterogeneous example

plot_pair <- function(evarname_str, pvarname_str, estimate_limits=c(0.01,.35)) {
  test_1 <- expos_wide |> filter(evarname == evarname_str, pvarname == pvarname_str) |> select(Begin.Year,  exposure_table_name, phenotype_table_name, e_variable_description, p_variable_description, estimate_adjusted, std.error_adjusted, p.value_adjusted) 
  exposure_name <- remove_units_from_string(test_1$e_variable_description[1])
  phenotype_name <- remove_units_from_string(test_1$p_variable_description[1])
  test_1 <- test_1 |> select(Begin.Year, estimate_adjusted, std.error_adjusted) |> rename(estimate=estimate_adjusted, std.error = std.error_adjusted, Survey=Begin.Year) |> mutate(i.squared = NA, i.squared_text='')
meta_test_1 <- adjusted_meta_2 |> filter(evarname == evarname_str, pvarname == pvarname_str) |> mutate(Survey = 'overall') |> select(Survey, estimate, std.error, i.squared) |> mutate(i.squared_text = sprintf("%i%%", round(i.squared)))
  test_1 <- test_1 |> rbind(meta_test_1) 
  test_1 <- test_1 |> mutate(point_shape = ifelse(Survey == 'overall', 23, 21)) |> mutate(point_shape = as.integer(point_shape)) 
  test_1 <- test_1 |> mutate(point_size = ifelse(Survey == 'overall', 7, 2)) 
  p <- ggplot(test_1, aes(Survey, estimate))
  p <- p + geom_point(aes(shape=point_shape, size=point_size, fill=point_shape)) + scale_shape_identity() + scale_size_identity()
  p <- p + geom_text(data=meta_test_1, aes(Survey, estimate, label=i.squared_text), size=3, color="black", nudge_x=.6) 
  p <- p + geom_errorbar(aes(ymin=estimate-1.96*std.error, ymax=estimate+1.96*std.error), width=.1) + scale_x_discrete(limits=rev)
  p <- p + scale_y_continuous(limits=estimate_limits)
  p <- p + coord_flip()  + theme_bw() + theme(legend.position = "none")
  p <- p + ggtitle(sprintf('scale(%s)-scale(log10(%s))', phenotype_name, exposure_name))+ theme(plot.title = element_text(size = 7))
}

## non-heterogeneous example
p1 <- plot_pair('URXP01', 'LBDNENO')
## heterogeneous example
p2 <- plot_pair('LBXGTC', 'BMXBMI')

p3 <- plot_pair('LBXBPB', 'BMXHT', c()) # 33% i2

p4 <- plot_pair('LBXCOT', 'BPXPLS', c()) # 33% i2



#expos_wide |> filter(evarname == 'LBXPFOS', pvarname == 'LBXSAL') 
plot_grid(p1, p2, p3, p4, ncol=2,labels = c('A', 'B', "C", "D"), label_size = 12)

# Examples for paper of top hits

rbind(
  adjusted_meta_2 |> filter(evarname == 'LBXPFHS', pvarname == 'LBXSAL') |> select(tau.squared),
  adjusted_meta_2 |> filter(evarname == 'LBXGTC', pvarname == 'LBXTC') |> select(tau.squared),
  adjusted_meta_2 |> filter(evarname == 'LBXBPB', pvarname == 'BMXWT') |> select(tau.squared),
  adjusted_meta_2 |> filter(evarname == 'LBXCOT', pvarname == 'LBDNENO') |> select(tau.squared)
)

## # A tibble: 4 × 1
##   tau.squared
##         <dbl>
## 1  0.0000737 
## 2  0         
## 3  0.0000634 
## 4  0.00000143

adjusted_meta_2 |> filter(evarname == 'LBXPFHS') |> group_by(sig_levels) |> summarize(sd_estimate=sd(estimate))

## # A tibble: 3 × 2
##   sig_levels   sd_estimate
##   <fct>              <dbl>
## 1 > BY & Bonf.      0.0282
## 2 BY<0.05           0.0682
## 3 Bonf.<0.05        0.107

adjusted_meta_2 |> filter(evarname == 'LBXGTC') |> group_by(sig_levels) |> summarize(sd_estimate=sd(estimate))

## # A tibble: 3 × 2
##   sig_levels   sd_estimate
##   <fct>              <dbl>
## 1 > BY & Bonf.      0.0269
## 2 BY<0.05           0.0587
## 3 Bonf.<0.05        0.143

adjusted_meta_2 |> filter(evarname == 'LBXBPB') |> group_by(sig_levels) |> summarize(sd_estimate=sd(estimate))

## # A tibble: 3 × 2
##   sig_levels   sd_estimate
##   <fct>              <dbl>
## 1 > BY & Bonf.      0.0349
## 2 BY<0.05           0.0566
## 3 Bonf.<0.05        0.0858

adjusted_meta_2 |> filter(evarname == 'LBXCOT') |> group_by(sig_levels) |> summarize(sd_estimate=sd(estimate))

## # A tibble: 3 × 2
##   sig_levels   sd_estimate
##   <fct>              <dbl>
## 1 > BY & Bonf.      0.0268
## 2 BY<0.05           0.0441
## 3 Bonf.<0.05        0.0814

p1 <- plot_pair('LBXPFHS', 'LBXSAL') + ylab("Association Estimate")
## heterogeneous example
p2 <- plot_pair('LBXGTC', 'LBXTC') + ylab("Association Estimate")

p3 <- plot_pair('LBXBPB', 'BMXWT', c(-.3, 0)) + ylab("Association Estimate")

p4 <- plot_pair('LBXCOT', 'LBDNENO')+ ylab("Association Estimate")

plot_grid(p1, p2, p3, p4, ncol=2,labels = c('A', 'B', "C", "D"), label_size = 12)

library(ggridges)
## show histogram of associations for all top findings for  LBXGTC, "LBXBPB", "LBXPFHS", "LBXCOT"
ep_candidates <-  tibble(evarname = c("LBXGTC", "LBXBPB", "LBXPFHS", "LBXCOT"), 
                         pvarname = c("LBXTC", "BMXWT", "LBXSAL", "LBDNENO"))
exposure_dist <- adjusted_meta_2 |> filter(evarname %in% ep_candidates$evarname) |> filter(sig_levels == "Bonf.<0.05") |> mutate(evardesc = remove_units_from_string(evardesc))

## collect Survey specific associations
survey_exposure_pts <- vector(mode = "list", length=nrow(ep_candidates))
for(rw_num in 1:nrow(ep_candidates)) {
  survey_exposure_pts[[rw_num]] <- expos_wide |> filter(evarname == ep_candidates$evarname[rw_num], pvarname == ep_candidates$pvarname[rw_num]) |> 
  select(Begin.Year,  evarname, pvarname, exposure_table_name, phenotype_table_name, e_variable_description, p_variable_description, estimate_adjusted, std.error_adjusted, p.value_adjusted) 
}

survey_exposure_pts <- survey_exposure_pts |> bind_rows() |> mutate(evardesc=remove_units_from_string(e_variable_description))

exposure_dist <- exposure_dist |> mutate(evarname = substr(evarname, 4, 8)) 
survey_exposure_pts <- survey_exposure_pts |> mutate(pvarname = substr(pvarname, 4, 8), evarname = substr(evarname, 4, 8))

p <- ggplot(exposure_dist, aes(y=evarname, x=estimate, fill=evarname))
p <- p + geom_density_ridges()+ scale_fill_jama(guide="none") #+ scale_colour_continuous(guide = "none")
p <- p + geom_point(data=survey_exposure_pts, aes(y=evarname, x=estimate_adjusted, color=factor(Begin.Year))) + scale_color_aaas(name = "")
pheno_het <- p + theme_bw()  + theme(legend.position='bottom') + ylab("Exposure Name") + xlab("Association Estimate")

survey_het <- plot_grid(p1, p2, p3, p4, ncol=2,labels = c('A', 'B', "C", "D"), label_size = 12)

p <- plot_grid(survey_het, pheno_het, ncol=2, labels=c("", "E"), label_size = 12)

## Picking joint bandwidth of 0.0395

p <- ggplot(adjusted_meta_2 , aes(estimate, estimate_unadjusted, color=sig_levels)) ## have to get this
p <- p + geom_point(shape='.') + scale_x_continuous(limits=c(-1, 1)) + scale_y_continuous(limits=c(-1, 1)) + scale_color_aaas()
p <- p + geom_abline()
p <- p + facet_grid(~sig_levels) + xlab("Adjusted model estimate [Exposure + Demographics]") + ylab("Unadjusted estimate [Exposure]")
p <- p + geom_smooth(method="lm")
p <- p + theme_bw() +theme(legend.position = "none")
p

## `geom_smooth()` using formula 'y ~ x'

## Warning: Removed 193 rows containing non-finite values (stat_smooth).

## Warning: Removed 193 rows containing missing values (geom_point).

# 
# 
# 
tidy(lm(estimate ~ estimate_unadjusted, adjusted_meta_2))

## # A tibble: 2 × 5
##   term                estimate std.error statistic   p.value
##   <chr>                  <dbl>     <dbl>     <dbl>     <dbl>
## 1 (Intercept)         -0.00471  0.000177     -26.6 1.27e-155
## 2 estimate_unadjusted  0.464    0.00141      329.  0

tidy(lm(estimate ~ estimate_unadjusted, adjusted_meta_2 |> filter(pvalue_bonferroni < .05)))

## # A tibble: 2 × 5
##   term                estimate std.error statistic  p.value
##   <chr>                  <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)          -0.0153  0.000877     -17.4 6.27e-67
## 2 estimate_unadjusted   0.631   0.00476      133.  0

tidy(lm(estimate ~ estimate_unadjusted, adjusted_meta_2 |> filter(sig_levels == '> BY & Bonf.')))

## # A tibble: 2 × 5
##   term                estimate std.error statistic  p.value
##   <chr>                  <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)         -0.00350  0.000180     -19.5 2.26e-84
## 2 estimate_unadjusted  0.421    0.00155      272.  0

tidy(lm(estimate ~ estimate_unadjusted, adjusted_meta_2 |> filter(sig_levels == 'BY<0.05')))

## # A tibble: 2 × 5
##   term                estimate std.error statistic  p.value
##   <chr>                  <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)         -0.00639  0.000632     -10.1 6.38e-24
## 2 estimate_unadjusted  0.485    0.00460      105.  0

# 
p <- ggplot(adjusted_meta_2 , aes(estimate_unadjusted, estimate, color=ecategory))
p <- p + geom_point(shape='.') + scale_x_continuous(limits=c(-1, 1)) + scale_y_continuous(limits=c(-1, 1)) + scale_color_aaas()
p <- p + geom_abline()
p <- p + facet_grid(~ecategory) + ylab("Adjusted model estimate [Exposure + Demographics]") + xlab("Unadjusted estimate [Exposure]")
p <- p + geom_smooth(method="lm")
p <- p + theme_bw() +theme(legend.position = "none")
p

## `geom_smooth()` using formula 'y ~ x'

## Warning: Removed 193 rows containing non-finite values (stat_smooth).
## Removed 193 rows containing missing values (geom_point).

Which e and P domains are most subject to demographic bias?

• nutrients and pollutants

bias_per_ecategory <- adjusted_meta_2 |> group_by(ecategory) |> summarize(
  mod=list(lm(estimate ~ estimate_unadjusted))) |> mutate(tidied=map(mod, tidy)) |> unnest(tidied)

bias_per_ecategory |> select(-mod) |> select(ecategory, term, estimate, p.value) |> filter(term == "estimate_unadjusted") |> select(-term)

## # A tibble: 9 × 3
##   ecategory         estimate  p.value
##   <chr>                <dbl>    <dbl>
## 1 allergy              0.753 2.65e-94
## 2 autoantibody         0.478 3.09e-64
## 3 housing              0.354 4.72e-81
## 4 income               0.725 0       
## 5 infection            0.557 0       
## 6 nutrients            0.287 0       
## 7 physical activity    0.415 0       
## 8 pollutant            0.386 0       
## 9 smoking              0.660 0

bias_per_pcategory <- adjusted_meta_2 |> group_by(pcategory) |> summarize(
  mod=list(lm(estimate ~ estimate_unadjusted))) |> mutate(tidied=map(mod, tidy)) |> unnest(tidied)

bias_per_pcategory |> select(-mod) |> select(pcategory, term, estimate, p.value) |> filter(term == "estimate_unadjusted") |> select(-term)

## # A tibble: 8 × 3
##   pcategory      estimate   p.value
##   <chr>             <dbl>     <dbl>
## 1 aging             0.300 9.28e- 16
## 2 anthropometric    0.347 0        
## 3 biochemistry      0.665 0        
## 4 blood             0.675 0        
## 5 blood pressure    0.514 0        
## 6 fitness           0.637 3.56e-159
## 7 lung              0.379 8.71e-193
## 8 microbiome        0.791 0

Multivariate R2 of the exposome

p <- ggplot(mvr2, aes(n_evars, mv*100))
p <- p + geom_point() + theme_bw()
p <- p + geom_text_repel(data=mvr2 |> filter(mv > .1),aes(n_evars, mv*100, label=substr(pvarname, 4, 10) ), max.overlaps = 20)
p <- p + xlab("Number of Exposome Variables in Model") + ylab("R-squared (%)")
p1 <- p

p <- ggplot(mvr2, aes(base_adj_rsq, mve_rsq))
p <- p + geom_point() + theme_bw()
p <- p + geom_text_repel(data=mvr2 |> filter(mv > .1),aes(base_adj_rsq, mve_rsq, label=substr(pvarname, 4, 10) ), max.overlaps = 20)
p <- p + xlab("R2 (Demographic Model)") + ylab("R2 (Full Model)")
p2 <- p

p3 <- ggplot(mvr2, aes(n_evars, mve_rsq))
p3 <- p3 + geom_point(aes(n_evars, base_adj_rsq))
p3 <- p3 + geom_segment(aes(x=n_evars, y=base_adj_rsq, xend=n_evars, yend=mve_rsq),arrow = arrow(length = unit(0.01, "npc")))
p3 <- p3 + geom_text_repel(data=mvr2 |> filter(mv > .1), aes(n_evars, mve_rsq, label=substr(pvarname, 4, 10) ), max.overlaps = 20)
p3 <- p3 + theme_bw()
p3

mvr2 |> summarize(mean_sample_size=mean(n), median_r2=median(mv), q25_r2=quantile(mv, probs=.25), q75_r2=quantile(mv, probs=.75))

## # A tibble: 1 × 4
##   mean_sample_size median_r2  q25_r2 q75_r2
##              <dbl>     <dbl>   <dbl>  <dbl>
## 1           16621.    0.0153 0.00458 0.0365

mvr2 |> summarize(n_min=min(n), n_max=max(n), n_evars_min=min(n_evars), n_evars_max=max(n_evars))

## # A tibble: 1 × 4
##   n_min n_max n_evars_min n_evars_max
##   <int> <int>       <int>       <int>
## 1  1602 68340           0          86

mvr2 |> arrange(-n_evars) |> filter(mv > .1) |> mutate(pvarname=substr(pvarname, 4, 8)) |> select(pvarname, pvardesc)  |> unite("index", pvarname, pvardesc, sep=":")  |> print(n=26)

## # A tibble: 16 × 1
##    index                                   
##    <chr>                                   
##  1 TRPF:Trunk Percent Fat                  
##  2 RBF:Folate, RBC (ng/mL RBC)             
##  3 SF1SI:5-Methyl-tetrahydrofolate (nmol/L)
##  4 FOT:Serum total folate (ng/mL)          
##  5 RFO:RBC folate (ng/mL)                  
##  6 SATA:Subcutaneous fat area              
##  7 SATM:Subcutaneous fat mass              
##  8 SATV:Subcutaneous fat volume            
##  9 TATA:Total abdominal fat area           
## 10 TATM:Total abdominal fat mass           
## 11 TATV:Total abdominal fat volume         
## 12 FOL:Folate, serum (ng/mL)               
## 13 STFAT:Subtotal Fat (g)                  
## 14 TOFAT:Total Fat (g)                     
## 15 TRFAT:Trunk Fat (grams)                 
## 16 THICR:Thigh Circumference (cm)

p <- plot_grid(p1, p2, ncol=2, labels=c("A", "B"))

Correlation of phenotypes in exposome space

library(corrr)
library(gplots)

## 
## Attaching package: 'gplots'

## The following object is masked from 'package:stats':
## 
##     lowess

to_array <- adjusted_meta_2 |> filter(expo_name == 'expo', vartype == 'continuous') |> select(evarname, pvarname, estimate, p.value) |> mutate(estimate= ifelse(p.value >= 0.05, 0, estimate)) |> mutate(estimate = ifelse(is.na(estimate), 0, estimate))  |> select(-p.value) |> pivot_wider(names_from = pvarname, values_from = estimate)

exposure_correlation <- to_array |> select(-evarname) |> correlate(diagonal = 1)

## Warning in stats::cor(x = x, y = y, use = use, method = method): the standard
## deviation is zero

## Correlation computed with
## • Method: 'pearson'
## • Missing treated using: 'pairwise.complete.obs'

m <- exposure_correlation |> as_matrix()
m[is.na(m)] <- 0

heatmapColors <- function(numColors=16) {
    c1 <- rainbow(numColors,v=seq(0.5,1,length=numColors),s=seq(1,0.3,length=numColors),start=4/6,end=4.0001/6);
    c2 <- rainbow(numColors,v=seq(0.5,1,length=numColors),s=seq(1,0.3,length=numColors),start=1/6,end=1.0001/6);
    c3 <- c(c1,rev(c2)); 
    return(c3)
}

heatmap.2(m, trace = 'none', na.rm = F, scale = 'none', symm = T, col=heatmapColors(5), margins=c(16,16), sepwidth=c(.1, .1), symbreaks=T)

Correlation of the exposome

to_array <- adjusted_meta_2 |> filter(expo_name == 'expo', vartype == 'continuous') |> select(evarname, pvarname, estimate, p.value) |> mutate(estimate= ifelse(p.value >= 0.05, 0, estimate)) |> mutate(estimate = ifelse(is.na(estimate), 0, estimate))  |> select(-p.value) |> pivot_wider(names_from = evarname, values_from = estimate)

phenome_correlation <- to_array |> select(-pvarname) |> correlate(diagonal = 1)

## Warning in stats::cor(x = x, y = y, use = use, method = method): the standard
## deviation is zero

## Correlation computed with
## • Method: 'pearson'
## • Missing treated using: 'pairwise.complete.obs'

m <- phenome_correlation |> as_matrix()
m[is.na(m)] <- 0

heatmap.2(m, trace = 'none', na.rm = F, scale = 'none', symm = T, col=heatmapColors(5), margins=c(16,16), sepwidth=c(.1, .1), symbreaks=T)